
function xOut=power_gen(x1,x2)
  n1 = size(x1);
  n2 = size(x2);
  x = getval(x1);
  y = getval(x2);
  dfdx1 = x.^(y-1).*y;
  dfdx2 = x.^y.*log(x);
  if(notisa(x1,'der1'))
    np = nindep(x2);
  else
    np = nindep(x1);
  end
  if(notisa(x1,'der1'))
    xOut.v = power_gen(x1,x2.v);
    xOut.d = elelmult_firstindx(dfdx2,x2.d);
  elseif(notisa(x2,'der1'))
    xOut.v = power_gen(x1.v,x2);
    xOut.d = elelmult_firstindx(dfdx1,x1.d);
  else  % both deriv1
    n1d = size(x1.d);
    np = n1d(end);
    xOut.v = power_gen(x1.v,x2.v);
    xOut.d = elelmult_firstindx(dfdx1,x1.d) + elelmult_firstindx(dfdx2,x2.d);
  end
  xOut=class(xOut,'der1');
